#!/usr/bin/perl

use strict;
use Getopt::Long;
use IO::File;

my $mode = "CC";
my $CCsamplefile = "/Net/spud/raid/CC-Cousin-Lines-SNP6/CClineIDs.csv";
my $inputdir = "/Net/spud/raid/CC-Cousin-Lines-SNP6/ANALYSIS-070809/genotype.10122009";
my $error = 0.001;
my $prefix="outbreds";
#my $MDAdir = "/Net/spud/raid/AFFY/MouseDiversityArray/";
my $snpmap = "/Net/spud/raid/CC-Cousin-Lines-SNP6/snpmap.txt";
my $selectfile = "outbred.ids.14122009.txt";

GetOptions( "dir=s"=>\$inputdir, "prefix=s"=>\$prefix, "select=s"=>\$selectfile );

my $selected;
my $family;
my $isfamily;
open(SELECT,$selectfile) || die "Could not open $selectfile\n";
while(<SELECT>) {
  chomp;
  my @a = split;
  $selected->{"$a[0].CEL"} = $a[1];
  $family->{"$a[0].CEL"} = $a[1];
  $isfamily->{$a[1]} = 1;
}

my $sampleinfo;
open(CC,$CCsamplefile)|| die "Could not open $CCsamplefile\n";
$_ = <CC>;
while(<CC>) {
  print;
  chomp;
  my ( $Position,$run,$CCline,$Generation,$Sex,$subline,$coatcolor,$run,$ID) = split(/,/);
  $sampleinfo->{$ID} = [$Position,$run,$CCline,$Generation,$Sex,$subline,$coatcolor,$run,$ID];
  print "$ID = [$Position,$run,$CCline,$Generation,$Sex,$subline,$coatcolor,$run,$ID]\n";
}

my ( $snpinfo, $genotype, $igenotype, $map, $jaxmap ) = ReadSNPMap($snpmap);

foreach my $chr ( sort keys %$map ) {
  my @sorted = sort { $a->[1] <=> $b->[1] } @{$map->{$chr}};
  $map->{$chr} = \@sorted;
}
my $genofile = "$inputdir/quant-norm.pm-only.brlmm-p.calls.txt";
open(GENO,$genofile) || die "Could not open $genofile\n";
print "reading $genofile\n";
my $n=0;
my @samples;
my $use;
my $geno;
my $nuse=0;

while (<GENO>) {
  chomp;
  if ( ! /^\#/ ) {
    if ( $n++ == 0 ) {
      my $dummy;
      ($dummy,@samples) = split;
      print "samples @samples\n";
      for (my $k=0;$k<=$#samples;$k++) {
	if ( $selected->{$samples[$k]} ) {
	  $use->[$k] = 1;
	  $nuse++;
	  print "$nuse use $k $samples[$k] $family->{$samples[$k]}\n";
	}	
      }
    } else {
      my ( $jax, @g ) = split;
      print "$n $jax\n" if ( !($n % 10000));
      for (my $k=0;$k<=$#g;$k++) {
	if ( $use->[$k] ){
	  $geno->{$samples[$k]}{$jax} = $genotype->{$jax}[$g[$k]+1];
	}
      }	
    }
  }
}

	
foreach my $chr ( sort keys %$map ) {
  my $handle;
  foreach my $family ( keys %$isfamily) {
    my $pedfile = "$chr.$family.ped";
    $handle->{$family} = new IO::File;
    open($handle->{$family}, ">$pedfile");
    print "writing $pedfile\n";
  }
  foreach my $sample ( sort keys %$geno ) {
    my $id = $sample;
    my $sex = "M";
    my $fam = $family->{$sample} ? $family->{$sample} : "CC";
    if ( $sampleinfo->{$id} ){
      $id = $sampleinfo->{$id}[2] . "-" . $sampleinfo->{$id}[5];
      $sex = $sampleinfo->{$id}[4];
      $fam = $sampleinfo->{$id}[2];
    }	
    $id =~ s/.CEL$//;
    my $txt = "$fam $id 0 0 M NA ";
    
    print {$handle->{$fam}} $txt;
    my $ns=0;
    my $na=0;
    foreach my $ref ( @{$map->{$chr}} ) {
      my $jax = $ref->[0];
      my $txt = $geno->{$sample}{$jax}? $geno->{$sample}{$jax} : " NA NA";
      print {$handle->{$fam}}  $txt;
      if ( $geno->{$sample}{$jax} ) {
	$ns++;
      }
      else {
	$na++;
      }	
    }	
    print {$handle->{$fam}}  "\n";
    my $sum = $ns+$na;
    print "chr $sample $ns $na $sum\n";
  }


  my $mapfile = "$chr.$prefix.map";
  open(MAP,">$mapfile") || die "Could not open $mapfile\n";
  print MAP "marker\tchromosome\tbp\n";



  foreach my $ref ( @{$map->{$chr}} ) {
    my $jax = $ref->[0];
    print MAP join( "\t", $jax, $chr, $ref->[1]) . "\n";
  }
  close(MAP);
}

sub ReadSNPMap {
  my ($snpmap) = @_;
  my $snpinfo;
  my $genotype;
  my $igenotype;
  my $map;
  my $jaxmap;

  open(SNPMAP,$snpmap) || die "Could not open $snpmap\n";
  print "reading $snpmap\n";
  my $nsnp=0;
  $_ = <SNPMAP>;
  while (<SNPMAP>) {
    chomp;
    my($jax,@data) = split(/,/);
    $snpinfo->{$jax} = \@data;
    $genotype->{$jax}[0] = "  NA NA";
    $genotype->{$jax}[1] = "  $data[3] $data[3]";
    $genotype->{$jax}[2] = "  $data[3] $data[4]";
    $genotype->{$jax}[3] = "  $data[4] $data[4]";
#    print "$jax $genotype->{$jax}[0] $genotype->{$jax}[1] $genotype->{$jax}[2] $genotype->{$jax}[3]\n";
    $igenotype->{$jax}[0] = undef;
    $igenotype->{$jax}[1] = $data[3].$data[3];
    $igenotype->{$jax}[3] = $data[4].$data[4];

    my $chr = $data[1];
    my $bp = $data[2];
    push @{$map->{$chr}}, [$jax, $bp];
    $jaxmap->{$jax} = [$chr, $bp];
    $nsnp++;
  }
  close(SNPMAP);
  print "read $nsnp snps\n";
  print "chrs " . join( " ", sort keys %$map) . "\n";
  return( $snpinfo, $genotype, $igenotype, $map, $jaxmap );

}

